function [mu] = MODIPcomp(lat, lon, stModip)
%
% Modified Dip Latitude (MODIP) Computation
%
%DESCRIPTION:
%This function computes the Modified Dip Latitude (MODIP, mu), 
%interpolating the array of MODIP grid.
%
%PROTOTYPE:
% [mu] = MODIPcomp(phi, lam, stModip)
%
%--------------------------------------------------------------------------
% INPUTS:
%   lat        [1x1]       Latitude                  [deg]
%   lon        [1x1]       Longitude                 [deg]
%   stModip    [39x39]     Array of MODIP Values     [-]
%--------------------------------------------------------------------------
% OUTPUTS:
%   mu         [1x1]       MODIP                     [deg]
%--------------------------------------------------------------------------
%
%NOTES:
% - Latitude should be between -90deg and +90deg
%
%CALLED FUNCTIONS:
% (none)
%
%UPDATES:
% (none)
%
%REFERENCES:
% [1] "Ionospheric Correction Algorithm for Galileo Single-Frequency Users"
%      - European GNSS (Galileo) Open Service
% [2] "Electron Density Models and Data for Transionospheric Radio
%      Propagation" - Report ITU-R P.2297-1 (05/2019)
%
%AUTHOR(s):
%Luigi De Maria, Matteo D'Addazio, 2022
%

%% Main Code

%Latitude Check
if lat >= +90       %Extreme Case 1
    mu = 90;
    
elseif lat <= -90   %Extreme Case 2
    mu = -90;
    
else                %Main Computations (-90deg < lat < +90deg)
    %First Index Determination
    l = round((lon + 180)/10) - 2;
    % l = round((lon + 180)/10) - 1;
    if     l < 0
        l = l + 36;
    elseif l > 36 %33
        l = l - 36;
    end

    
    %Initial Position Determination  %OBS This code has been revised w/ Annex F [1]
    a = (lat + 90)/5 + 1;
    
    %Second Index Determination
    x = a - round(a);
    % i = int32(a) - 2 - 1; %(the extra "-1" is to account for Python indexing from "-1")
    i = round(a) - 2;

    %Subgrid Extraction from MODIP Array
    z = zeros(4,4);
    % shift = 4;
    for k = 1 : 4
        for j = 1 : 4
            %if i+j>0 && l+k>0
                z(j,k) = stModip(i+j,l+k);
            % end
            % if i+j>0 && l+k<=0
            %     z(j,k) = stModip(i+j,end+(l+k+2));
            % end
            % if i+j<=0 && l+k >0
            %     z(j,k) = stModip(end+(i+j+2), l+k);
            % end
            % if i+j<=0 && l+k <=0
            %     z(j,k) = stModip(end+(i+j), end+(l+k+2));
            % end
        end
    end
    
    %Coordinates Interpolation
    zk = zeros(4,1);
    for k = 1 : 4
        zk(k) = interp3(z(1,k),z(2,k),z(3,k),z(4,k),x); 
    end
    
    %New Position Computation
    b = (lon + 180)/10;
    y = b - round(b);
    
    %MODIP [deg]
    mu = interp3(zk(1),zk(2),zk(3),zk(4),y);
end

end